# ------------------------------------------------------------------------------
# Script: Twisting-Induced Instabilities in Double-Helix Chiral Rods
# Author: Giada Risso
# Affiliation: SEAS, Harvard University
# Date: 2025-06-03
# Description: Automates parametric generation, meshing, boundary conditions,
#              and simulation setup for soft cylindrical helical structures.
# ------------------------------------------------------------------------------

# --- Import ABAQUS Modules ---
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *

import os
import sys
import time
import csv
import math
import subprocess
import regionToolset
import displayGroupMdbToolset as dgm
import xyPlot
import displayGroupOdbToolset as dgo

### PARAMETERS ###

projectPath = r'C:\Users\gir796\OneDrive - Harvard University\Documents\Harvard\Instabilities of Soft Cylinders\DATA\videos\Final Code' #Put the absolute path of the folder of the project

number_of_simulation = 1

# Geometric parameters (lenghts in millimeters)

rint_table = [4.0] * number_of_simulation  #       (in millimeters)
rext_table = [10.0] * number_of_simulation  #      (in millimeters)
p_table = [10.0] * number_of_simulation #          (in millimeters/round)
alpha_deg_table =  [45] * number_of_simulation #  (in degrees)
l_table = [110]* number_of_simulation  #          (in millimeters)
tc_table = [5]* number_of_simulation  #           (in millimeters)

#

size_factor_table = [1.0] * number_of_simulation # To proportionnaly dowmsize or upsize the geometry

contact_table = [' '] * number_of_simulation #'Hard' or 'Exp' otherwise put ' '

damping_factor_table = [1e-08] * number_of_simulation #Positive float or None if no damping stabilization required

# BC parameters

stretching_table = [0.1] * number_of_simulation # l / l0 - 1
twisting_table = [6.28]  #in radians, with a + for closing and - for opening

#

mesh_size_table = [1.0] * number_of_simulation # (in millimeters)

cores_table = [6] * number_of_simulation # Even

data_file_name = 'SingleRun_' # Name of the record of the simulations

###
session.journalOptions.setValues(replayGeometry=COORDINATE, recoverGeometry=COORDINATE) # Allow to catch object through coordinates

setPath = projectPath
os.chdir(setPath)

sim_time_table = []
job_name_table = []

dataPath = os.path.join(projectPath, data_file_name + '.csv')


with open(dataPath, 'w') as f:
    writer = csv.writer(f, delimiter=';', lineterminator='\n')
    writer.writerow(['Index'] + [i for i in range(number_of_simulation)])
    writer.writerow(['Rint (mm)'] + rint_table)
    writer.writerow(['Rext (mm)'] + rext_table)
    writer.writerow(['Pitch (mm/round)'] + p_table)
    writer.writerow(['Alpha (deg)'] + alpha_deg_table)
    writer.writerow(['Lenght (mm)'] + l_table)
    writer.writerow(['Cap thickness (mm)'] + tc_table)
    writer.writerow(['Stretching'] + stretching_table)
    writer.writerow(['Twisting (rad)'] + twisting_table)
    writer.writerow(['Mesh size (mm)'] + mesh_size_table)
    writer.writerow(['Number of cores'] + cores_table)

for n in range(number_of_simulation):

    Mdb()
    model = mdb.models['Model-1']

    ### PARAMETERS ###

    size_factor = size_factor_table[n]

    rint = rint_table[n] * size_factor
    rext = rext_table[n] * size_factor
    p = p_table[n]
    alpha_deg = alpha_deg_table[n]
    l = l_table[n] * size_factor
    tc = tc_table[n] * size_factor

    contact = contact_table[n]
    damping_factor = damping_factor_table[n]

    stretching = stretching_table[n]
    twisting = twisting_table[n]

    mesh_size = mesh_size_table[n] * size_factor

    cores = cores_table[n]

    ###

    alpha = alpha_deg * math.pi / 180.0 # Convert degrees into radians

    ### PART ###

    # Helix
    model.ConstrainedSketch(name='__profile__', sheetSize=100.0)
    model.sketches['__profile__'].Spot(point=(0.0, 0.0))
    model.sketches['__profile__'].ArcByCenterEnds(center=(0.0, 0.0)
        , direction=CLOCKWISE, point1=(-rint * math.cos(alpha / 2), rint * math.sin(alpha / 2)),
        point2=(rint * math.cos(alpha / 2), rint * math.sin(alpha / 2)))
    model.sketches['__profile__'].ArcByCenterEnds(center=(0.0, 0.0)
        , direction=COUNTERCLOCKWISE, point1=(-rint * math.cos(alpha / 2), -rint * math.sin(alpha / 2)),
        point2=(rint * math.cos(alpha / 2), -rint * math.sin(alpha / 2)))
    model.sketches['__profile__'].ArcByCenterEnds(center=(0.0, 0.0)
        , direction=COUNTERCLOCKWISE, point1=(-rext * math.cos(alpha / 2), rext * math.sin(alpha / 2)),
        point2=(-rext * math.cos(alpha / 2), -rext * math.sin(alpha / 2)))
    model.sketches['__profile__'].ArcByCenterEnds(center=(0.0, 0.0)
        , direction=CLOCKWISE, point1=(rext * math.cos(alpha / 2), rext * math.sin(alpha / 2)),
        point2=(rext * math.cos(alpha / 2), -rext * math.sin(alpha / 2)))
    model.sketches['__profile__'].Line(point1=(-rint * math.cos(alpha / 2), rint * math.sin(alpha / 2)),
        point2=(-rext * math.cos(alpha / 2), rext * math.sin(alpha / 2)))
    model.sketches['__profile__'].Line(point1=(-rint * math.cos(alpha / 2), -rint * math.sin(alpha / 2)),
        point2=(-rext * math.cos(alpha / 2), -rext * math.sin(alpha / 2)))
    model.sketches['__profile__'].Line(point1=(rint * math.cos(alpha / 2), -rint * math.sin(alpha / 2)),
        point2=(rext * math.cos(alpha / 2), -rext * math.sin(alpha / 2)))
    model.sketches['__profile__'].Line(point1=(rint * math.cos(alpha / 2), rint * math.sin(alpha / 2)),
        point2=(rext * math.cos(alpha / 2), rext * math.sin(alpha / 2)))
    model.Part(dimensionality=THREE_D, name='Helix', type=
        DEFORMABLE_BODY)
    model.parts['Helix'].BaseSolidExtrude(depth=l - 2 * tc, pitch=p,
        sketch=model.sketches['__profile__'])
    del model.sketches['__profile__']

    # Cap 1
    model.ConstrainedSketch(name='__profile__', sheetSize=100.0)
    model.sketches['__profile__'].sketchOptions.setValues(
        decimalPlaces=3)
    model.sketches['__profile__'].CircleByCenterPerimeter(center=(
        0.0, 0.0), point1=(rext, 0.0))
    model.Part(dimensionality=THREE_D, name='Cap-1', type=
        DEFORMABLE_BODY)
    model.parts['Cap-1'].BaseSolidExtrude(depth=tc, sketch=
        model.sketches['__profile__'])
    del model.sketches['__profile__']

    # Cap 2
    model.Part(name='Cap-2', objectToCopy=
        model.parts['Cap-1'])
    model.parts['Cap-1'].setValues(space=THREE_D, type=
        DEFORMABLE_BODY)

    ### INSTANCES ###

    model.rootAssembly.DatumCsysByDefault(CARTESIAN)
    model.rootAssembly.Instance(dependent=ON, name='Helix-1', part=
        model.parts['Helix'])
    model.rootAssembly.Instance(dependent=ON, name='Cap-1-1', part=
        model.parts['Cap-1'])
    model.rootAssembly.Instance(dependent=ON, name='Cap-2-1', part=
        model.parts['Cap-2'])

    model.rootAssembly.translate(instanceList=('Helix-1', ),
        vector=(0.0, 0.0, tc))
    model.rootAssembly.translate(instanceList=('Cap-2-1', ),
        vector=(0.0, 0.0, l - tc))

    ### MERGE ###

    model.rootAssembly.InstanceFromBooleanMerge(domain=GEOMETRY,
        instances=(model.rootAssembly.instances['Cap-1-1'],
        model.rootAssembly.instances['Cap-2-1'],
        model.rootAssembly.instances['Helix-1']),
        keepIntersections=ON, name='Coil', originalInstances=SUPPRESS)

     ### MATERIAL ###

    model.Material(name='Rubber')
    model.materials['Rubber'].Elastic(table=((1.2, 0.495), ))

    ### SECTION ###

    model.HomogeneousSolidSection(material='Rubber', name=
        'Section-1', thickness=None)
    model.parts['Coil'].Set(cells=
        model.parts['Coil'].cells.getSequenceFromMask(('[#1f ]', ),
        ), name='Set-1')
    model.parts['Coil'].SectionAssignment(offset=0.0, offsetField=
        '', offsetType=MIDDLE_SURFACE, region=
        model.parts['Coil'].sets['Set-1'], sectionName='Section-1',
        thicknessAssignment=FROM_SECTION)
    
    ### MESH ###

    #Controls
    model.parts['Coil'].setMeshControls(elemShape=TET, regions=
        model.parts['Coil'].cells.getSequenceFromMask(('[#1f ]', ),
        ), technique=FREE)

    #Element type
    model.parts['Coil'].setElementType(elemTypes=(ElemType(
        elemCode=C3D8R, elemLibrary=STANDARD), ElemType(elemCode=C3D6,
        elemLibrary=STANDARD), ElemType(elemCode=C3D4H, elemLibrary=STANDARD,
        secondOrderAccuracy=OFF, distortionControl=DEFAULT)), regions=(
        model.parts['Coil'].cells.getSequenceFromMask(('[#1f ]', ),
        ), ))

    #Seeds
    model.parts['Coil'].seedPart(deviationFactor=0.5,
        minSizeFactor=0.1, size=mesh_size)

    model.parts['Coil'].generateMesh()

    ### STEPS ###

    #RP and constraints
    model.rootAssembly.ReferencePoint(point=
        model.rootAssembly.instances['Coil-1'].InterestingPoint(
        model.rootAssembly.instances['Coil-1'].edges.findAt((0.0,
        rext, 0.0), ), CENTER))
    model.rootAssembly.ReferencePoint(point=
        model.rootAssembly.instances['Coil-1'].InterestingPoint(
        model.rootAssembly.instances['Coil-1'].edges.findAt((0.0,
        -rext, l), ), CENTER))
    model.rootAssembly.Set(faces=
        model.rootAssembly.instances['Coil-1'].faces.findAt(((
        0.1030057 * rext, -0.9800334 * rext, 0.0), )), name='t_Set-1')
    model.RigidBody(name='Constraint-1', refPointRegion=Region(
        referencePoints=(model.rootAssembly.referencePoints[10], ))
        , tieRegion=model.rootAssembly.sets['t_Set-1'])
    model.rootAssembly.Set(faces=
        model.rootAssembly.instances['Coil-1'].faces.findAt(((
        0.2357023 * rext, 0.8516219 * rext, l), )), name='t_Set-2')
    model.RigidBody(name='Constraint-2', refPointRegion=Region(
        referencePoints=(model.rootAssembly.referencePoints[11], ))
        , tieRegion=model.rootAssembly.sets['t_Set-2'])

    e1 = 0.4 * mesh_size
    e2 = 0.75 * mesh_size
    central_line = model.parts['Coil'].nodes.getByBoundingBox(-e1, -e1, tc, e1, e1, l - tc)
    ridge_line = model.parts['Coil'].nodes.getByBoundingBox(rext - e2, -e2, tc, rext + e2, e2, l - tc)
    model.parts['Coil'].Set(nodes=central_line, name='Central-line')
    model.parts['Coil'].Set(nodes=ridge_line, name='Ridge-line')
   
    #Stretching
    model.StaticStep(initialInc=0.1, maxInc=0.1, maxNumInc=10000, 
        minInc=1e-12, name='Stretching', nlgeom=ON, previous='Initial')
    if (type(damping_factor) == float):
        if damping_factor > 0:
            model.steps['Stretching'].setValues(adaptiveDampingRatio=None,
                continueDampingFactors=False, stabilizationMagnitude=damping_factor,
                stabilizationMethod=DAMPING_FACTOR)
    model.EncastreBC(createStepName='Stretching', localCsys=None,
        name='Encastre', region=Region(referencePoints=(
        model.rootAssembly.referencePoints[10], )))
    model.DisplacementBC(amplitude=UNSET, createStepName=
        'Stretching', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=
        None, name='Stretching', region=Region(referencePoints=(
        model.rootAssembly.referencePoints[11], )), u1=0.0, u2=0.0,
        u3=stretching * l, ur1=0.0, ur2=0.0, ur3=0.0)

    #Twisting
    model.StaticStep(initialInc=0.1, maxInc=0.1, maxNumInc=10000,
        minInc=1e-12, name='Twisting', previous='Stretching')
    if (type(damping_factor) == float):
        if damping_factor > 0:
            model.steps['Twisting'].setValues(adaptiveDampingRatio=None,
                continueDampingFactors=False, stabilizationMagnitude=damping_factor,
                stabilizationMethod=DAMPING_FACTOR)
    model.boundaryConditions['Stretching'].deactivate('Twisting')
    model.DisplacementBC(amplitude=UNSET, createStepName='Twisting'
        , distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None, name=
        'Twisting', region=Region(referencePoints=(
        model.rootAssembly.referencePoints[11], )), u1=0.0, u2=0.0,
        u3=stretching * l, ur1=0.0, ur2=0.0, ur3=twisting)

    ### CONTACT ###

    #Create the surface (here whole)
    a = model.rootAssembly
    s1 = a.instances['Coil-1'].faces
    side1Faces1 = s1.getSequenceFromMask(mask=('[#2fdff ]', ), )
    region=a.Surface(side1Faces=side1Faces1, name='Surf-1')

    #Hard contact
    if contact == 'Hard':
        model.ContactProperty('Hard')
        model.interactionProperties['Hard'].TangentialBehavior(
            formulation=FRICTIONLESS)
        model.interactionProperties['Hard'].NormalBehavior(
            pressureOverclosure=HARD, allowSeparation=ON, 
            constraintEnforcementMethod=DEFAULT)

        model.SelfContactStd(name='Contact-Hard', createStepName='Initial', 
            surface=region, interactionProperty='Hard', thickness=ON)

    #Exponential contact
    if contact == 'Exp':
        model.ContactProperty('Exp')
        model.interactionProperties['Exp'].TangentialBehavior(
            formulation=FRICTIONLESS)
        model.interactionProperties['Exp'].NormalBehavior(
            pressureOverclosure=EXPONENTIAL, table=((2.0, 0.0), (0.0, 0.0022)), 
            maxStiffness=None, constraintEnforcementMethod=DEFAULT)

        model.SelfContactStd(name='Contact-Exp', createStepName='Initial', 
        surface=region, interactionProperty='Exp', thickness=ON)

    ### OUTPUTS REQUEST ###

    model.fieldOutputRequests['F-Output-1'].setValuesInStep(stepName='Stretching', variables=('S', 'LE', 'U', 'RF', 'COORD'))

    ### JOB ###

    def twist_direction(x):
        if x >= 0:
            return 'closing'
        else:
            return 'opening'

    job_name =  ('Job_' + data_file_name + '_rint=' + str(rint) + '_rext=' + str(rext) + '_p=' + str(round(p, 1)) + '_a=' + str(alpha_deg) + '_s=' + str(stretching) + '_' + twist_direction(twisting)).replace('.', ',') # To choose (be careful of the lenght of the name)

    mdb.Job(atTime=None, contactPrint=OFF, description='', echoPrint=OFF,
        explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF,
        memory=90, memoryUnits=PERCENTAGE, model='Model-1', modelPrint=OFF,
        multiprocessingMode=DEFAULT, name=job_name, nodalOutputPrecision=SINGLE,
        numCpus=cores, numDomains=cores, numGPUs=0, queue=None, resultsFormat=ODB, scratch=
        '', type=ANALYSIS, userSubroutine='', waitHours=0, waitMinutes=0)
    start_time = time.time()
    # uncomment this line to run the job automatically
    # mdb.jobs[job_name].submit(consistencyChecking=OFF)
    mdb.jobs[job_name].waitForCompletion()
    end_time = time.time()
    sim_time = end_time - start_time
    sim_time_table.append(sim_time)
    job_name_table.append(job_name)